Geospatial Carpentry | Session 4
Clémentine Cottineau
7/18/2022
Let’s import important libraries
library(dplyr)
library(sf)
library(leaflet)
library(ggplot2)
library(tibble)
library(rgdal)
library(proj4)
library(remotes)
remotes::install_github('ropensci/osmdata')
library(osmdata)
library(osmextract)
# proj_db <- system.file("proj/proj.db", package = "sf")Import vector data from Open Street Map
# from : https://cran.r-project.org/web/packages/osmdata/vignettes/osmdata.html
bb <- getbb('rotterdam, NL', format_out = 'sf_polygon')
x <- opq(bbox = bb) %>%
add_osm_feature(key = 'building') %>%
osmdata_sf ()
str(x$osm_polygons)## Classes 'sf' and 'data.frame': 253568 obs. of 274 variables:
## $ osm_id : chr "23270254" "26545811" "34294334" "34294530" ...
## $ name : chr NA "Montevideo" "Lotus Serre" "Ingang" ...
## $ X3dr.height1 : chr NA NA NA NA ...
## $ X3dr.length1 : chr NA NA NA NA ...
## $ X3dr.type : chr NA NA NA NA ...
## $ X3dshapes.note : chr NA NA NA NA ...
## $ access : chr NA NA "customers" NA ...
## $ addr.city : chr NA NA NA NA ...
## $ addr.country : chr NA NA NA NA ...
## $ addr.full : chr NA NA NA NA ...
## $ addr.housename : chr NA NA NA NA ...
## $ addr.housenumber : chr NA NA NA NA ...
## $ addr.postcode : chr NA NA NA NA ...
## $ addr.street : chr NA NA NA NA ...
## $ addr.unit : chr NA NA NA NA ...
## $ aeroway : chr NA NA NA NA ...
## $ air_conditioning : chr NA NA NA NA ...
## $ alt_name : chr NA NA NA NA ...
## $ amenity : chr NA NA "restaurant" NA ...
## $ architect : chr NA NA NA NA ...
## $ area : chr NA NA NA NA ...
## $ artist_name : chr NA NA NA NA ...
## $ artwork_type : chr NA NA NA NA ...
## $ atm : chr NA NA NA NA ...
## $ bar : chr NA NA NA NA ...
## $ bicycle : chr NA NA NA NA ...
## $ bicycle_parking : chr NA NA NA NA ...
## $ brand : chr NA NA NA NA ...
## $ brand.wikidata : chr NA NA NA NA ...
## $ brand.wikipedia : chr NA NA NA NA ...
## $ bridge.support : chr NA NA NA NA ...
## $ building : chr "transportation" "apartments" "yes" "yes" ...
## $ building.architecture : chr NA NA NA NA ...
## $ building.colour : chr NA NA NA NA ...
## $ building.flats : chr NA NA NA NA ...
## $ building.levels : chr "2" "43" NA NA ...
## $ building.levels.underground : chr NA NA NA NA ...
## $ building.material : chr NA NA NA NA ...
## $ building.min_level : chr NA NA NA NA ...
## $ building.part : chr NA NA NA NA ...
## $ building.use : chr NA NA NA NA ...
## $ capacity : chr NA NA NA NA ...
## $ changing_table : chr NA NA NA NA ...
## $ check_date : chr NA NA NA NA ...
## $ check_date.opening_hours : chr NA NA NA NA ...
## $ cinema.3D : chr NA NA NA NA ...
## $ climbing.length.max : chr NA NA NA NA ...
## $ climbing.sport : chr NA NA NA NA ...
## $ clothes : chr NA NA NA NA ...
## $ club : chr NA NA NA NA ...
## $ community_centre : chr NA NA NA NA ...
## $ company : chr NA NA NA NA ...
## $ construction : chr NA NA NA NA ...
## $ contact.website : chr NA NA NA NA ...
## $ content : chr NA NA NA NA ...
## $ covered : chr NA NA NA NA ...
## $ craft : chr NA NA NA NA ...
## $ cuisine : chr NA NA NA NA ...
## $ delivery : chr NA NA NA NA ...
## $ demolished : chr NA NA NA NA ...
## $ denomination : chr NA NA NA NA ...
## $ denomination.nl : chr NA NA NA NA ...
## $ description : chr NA NA NA NA ...
## $ description.de : chr NA NA NA NA ...
## $ designation : chr NA NA NA NA ...
## $ dhm_id : chr NA NA NA NA ...
## $ disused : chr NA NA NA NA ...
## $ dock : chr NA NA NA NA ...
## $ drinking_water : chr NA NA NA NA ...
## $ email : chr NA NA NA NA ...
## $ emergency : chr NA NA NA NA ...
## $ facebook : chr NA NA NA NA ...
## $ fax : chr NA NA NA NA ...
## $ fee : chr NA NA NA NA ...
## $ ferry : chr NA NA NA NA ...
## $ fixme : chr NA NA NA NA ...
## $ floating : chr NA NA NA NA ...
## $ food : chr NA NA NA NA ...
## $ foot : chr NA NA NA NA ...
## $ frequency : chr NA NA NA NA ...
## $ fuel.GTL_diesel : chr NA NA NA NA ...
## $ fuel.diesel : chr NA NA NA NA ...
## $ fuel.lpg : chr NA NA NA NA ...
## $ fuel.octane_95 : chr NA NA NA NA ...
## $ fuel.octane_98 : chr NA NA NA NA ...
## $ generator.source : chr NA NA NA NA ...
## $ government : chr NA NA NA NA ...
## $ healthcare : chr NA NA NA NA ...
## $ healthcare.speciality : chr NA NA NA NA ...
## $ height : chr NA "152.3" NA NA ...
## $ heritage : chr NA NA NA NA ...
## $ heritage.operator : chr NA NA NA NA ...
## $ historic : chr NA NA NA NA ...
## $ image : chr NA "https://upload.wikimedia.org/wikipedia/commons/f/f2/Montevideo%28Rotterdam%29.jpg" NA NA ...
## $ image.0 : chr NA NA NA NA ...
## $ image.1 : chr NA NA NA NA ...
## $ image.2 : chr NA NA NA NA ...
## $ image.3 : chr NA NA NA NA ...
## $ indoor : chr NA NA NA NA ...
## [list output truncated]
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:273] "osm_id" "name" "X3dr.height1" "X3dr.length1" ...Map building age focusing on post-war buildings
buildings <- x$osm_polygons %>% st_transform(.,crs=28992)
buildings$build_date <- ifelse(buildings$start_date <1945, 1945,buildings$start_date)
ggplot(data = buildings) +
geom_sf(aes(fill = as.numeric(build_date), colour=as.numeric(build_date))) +
scale_fill_viridis_c(option = "viridis")+
scale_colour_viridis_c(option = "viridis")
Let’s focus on really old building and imagine we’re in charge of conservation. We want to know how much of the city would be affected by a non-construction zone of 500m around pre-1800 buildings.
Selecting old buildings
old <- 1800 # year prior to which you consider a building old
distance <- 500 # in meters
old_buildings <- filter(buildings, start_date <= old)
ggplot(data = old_buildings) + geom_sf()
Creating a buffer around a polygon
buffer_old_buildings <- st_buffer(old_buildings, dist = distance)
ggplot(data = buffer_old_buildings) + geom_sf()
Fusing overlapping buffers into a single polygon
single_old_buffer <- st_union(buffer_old_buildings) %>%
st_cast(to = "POLYGON") %>% st_as_sf() %>%
add_column("ID"=as.factor(1:nrow(.))) %>% st_transform(.,crs=4326) Representing old buildings by their centroid
sf::sf_use_s2(FALSE)
centroids_old <- st_centroid(old_buildings) %>% st_transform(.,crs=4326)
ggplot() +
geom_sf(data = single_old_buffer, aes(fill=ID)) +
geom_sf(data = centroids_old)
Counting the number of centroid per buffer with st_intersection (/! ≠st_intersect!)
centroids_buffers <- st_intersection(centroids_old,single_old_buffer) %>%
add_column(n=1)
centroid_by_buffer <- centroids_buffers %>%
group_by(ID) %>%
summarise(
n = sum(n)
)
single_buffer <- single_old_buffer %>% add_column(n_buildings = centroid_by_buffer$n)
ggplot() +
geom_sf(data = single_buffer, aes(fill=n_buildings)) +
scale_fill_viridis_c(alpha = 0.8,
begin = 0.6,
end = 1,
direction = -1,
option = "B")
Final visualisation
ggplot() +
geom_sf(data = buildings) +
geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
scale_fill_viridis_c(alpha = 0.6,
begin = 0.6,
end = 1,
direction = -1,
option = "B") 
old <- 1939 # year prior to which you consider a building old
distance <- 100 # in metersNow let’s look at the same analysis but for pre-war buildings and 100m buffers
old <- 1939
distance <- 100
#select
old_buildings <- filter(buildings, start_date <= old)
#buffer
buffer_old_buildings <- st_buffer(old_buildings, dist = distance)
#union
single_old_buffer <- st_union(buffer_old_buildings) %>%
st_cast(to = "POLYGON") %>% st_as_sf() %>%
add_column("ID"=1:nrow(.)) %>% st_transform(.,crs=4326)
#centroids
centroids_old <- st_centroid(old_buildings) %>% st_transform(.,crs=4326)
#intersection
centroids_buffers <- st_intersection(centroids_old,single_old_buffer) %>%
add_column(n=1)
centroid_by_buffer <- centroids_buffers %>%
group_by(ID) %>%
summarise(
n = sum(n)
)
single_buffer <- single_old_buffer %>% add_column(n_buildings = centroid_by_buffer$n)
ggplot() +
geom_sf(data = buildings) +
geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
scale_fill_viridis_c(alpha = 0.6,
begin = 0.6,
end = 1,
direction = -1,
option = "B") 
Problem: there are many pre-war buildings and the buffers are large so the number of old buildings is not very meaningful. Let’s compute the density of old buildings per buffer zone.
Visualisation
single_buffer$area <- st_area(single_buffer, ) %>% units::set_units(., km^2)
single_buffer$old_buildings_per_km2 <- as.numeric(single_buffer$n_buildings / single_buffer$area)
ggplot() +
geom_sf(data = buildings) +
geom_sf(data = single_buffer, aes(fill=old_buildings_per_km2), colour = NA) +
scale_fill_viridis_c(alpha = 0.6,
begin = 0.6,
end = 1,
direction = -1,
option = "B") 
EXERCISE!
Get OSM limits of Rotterdam
rdam <- getbb(
place_name = "Rotterdam",
format_out = 'sf_polygon',
silent = FALSE
) ## [1] "https://nominatim.openstreetmap.org/search?format=json&q=Rotterdam&featuretype=settlement&polygon_text=1&limit=10"Import and reproject CBS 500mx500m grid cells for the entire NL
# downloaded from: https://www.cbs.nl/-/media/cbs/dossiers/nederland-regionaal/vierkanten/500/2022-cbs_vk500_2019_vol.zip
cells_NL <- st_read("2019-cbs_vk500_2015_v2/CBS_VK500_2015_v2.shp") %>%
mutate_if(is.numeric, funs(ifelse(.==-99997, NA, .)))## Reading layer `CBS_VK500_2015_v2' from data source
## `/Users/ccottineau/Documents/Teaching/DataCarpentryWorkshop/GeoSpatial/2019-cbs_vk500_2015_v2/CBS_VK500_2015_v2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 151108 features and 131 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 13000 ymin: 306500 xmax: 278500 ymax: 619500
## Projected CRS: Amersfoort / RD Newcells_NL <- st_transform(cells_NL,crs=4326) INTERSECT (≠ st_intersection) to keep only cells for rotterdam
intersect_cells_rdam <- st_intersects(cells_NL, rdam, sparse=F)
intersect_cells_rdam_sf <- cbind(intersect_cells_rdam,cells_NL)
cells_rdam <- intersect_cells_rdam_sf %>%
filter(X1 == T )
plot(cells_rdam)